######
######
###### "Purging militaries : Introducing the Military Purges in Dictatorships (MPD) dataset 
###### 
###### JPR Replication file
######
###### 

## Jun Koga Sudduth 
### May 2020

library(tidyverse)
library(lubridate) # time date
library(foreign)
#install.packages("readstata13")
library(readstata13) 
#install.packages("feather")
library(feather)
library(broom)## model building
#install.packages("margins")
library(margins)
#install.packages("GLMMadaptive")
library(GLMMadaptive)
##install.packages("rms")
library(rms) 
library(gridExtra)

 

####
#####
#### Basic Model 
####
#####

setwd("C:/Users/Dell User/Dropbox/Military Purge Dataset Project")
dat <- read_csv("Purgedata_version2_2019/FINAL/analysis_tscs_v1.csv" ) 
## A tibble: 4,026 x 26
 

## Non-democracy  
dat <- dat %>%
  filter(gwf_regime!="democracy" )

dat$gwf_regime <- ordered(dat$gwf_regime, levels=c(  "personal","military", "party", "monarchy", "foreign/warload"))

 



######## 
######## Figure 1 (in Main-Text)
######## 

##### Trend of purges by region and year 

## country-year format 
## How many military purge cases (i.e. purge-year) per year 

p1 <- dat %>% filter(gwf_regime!="democracy" & gwf_regime!="foreign/warload")%>%  
  group_by(ccode, year) %>% # country-year format
  summarise(
    region2 = first(region2), 
    sum_purge = sum(purge_v1) 
  ) %>% ungroup() %>% 
  mutate( 
    purge = ifelse(sum_purge>0, 1, 0), 
    region2=ifelse(region2=="Eastern/Western Europe", "Eastern/Western \n Europe", 
                   ifelse(region2=="Latin America", "Latin \n America", 
                          ifelse(region2=="North Africa and the Middle East", "North Africa and \n Middle East", 
                                 ifelse(region2=="Sub-Saharan Africa", "Sub-Saharan \n Africa", region2))))
  ) %>% 
  filter(purge==1)%>%     # group_by(region2, year) %>% summarise(count=n()) %>% print(n=Inf)
  ggplot()+
  geom_bar(aes(year, fill=region2), width=0.4)+     ## the height of the bar proportional to the number of cases in each group
  facet_wrap(~ region2,  nrow = 1)   +
  ylab("Number of purges")+
  xlab("Year") +
  theme(legend.position="none") 



## Probability of Purges per Country-year 

p2 <- dat %>% filter(gwf_regime!="democracy" & gwf_regime!="foreign/warload")%>%  
  group_by(ccode, year) %>% 
  summarise(
    region2 = first(region2), 
    sum_purge = sum(purge_v1) , 
    leader= first(leader)
  ) %>% ungroup() %>% 
  mutate(   
    purge = ifelse(sum_purge>0, 1, 0) , 
    region2=ifelse(region2=="Eastern/Western Europe", "Eastern/Western \n Europe", 
                   ifelse(region2=="Latin America", "Latin \n America", 
                          ifelse(region2=="North Africa and the Middle East", "North Africa and \n Middle East", 
                                 ifelse(region2=="Sub-Saharan Africa", "Sub-Saharan \n Africa", region2))))  
     
  ) %>% 
  group_by(region2, year) %>% 
  summarise( 
    a_country = length( unique(ccode)),
    npurge = sum(purge, na.rm=TRUE), 
    purge_freq =npurge/ a_country, 
    fccode=first(ccode), 
    lccode=last(ccode), 
    fleader=first(leader), 
    lleader=last(leader) 
    
  )  %>%     
  ggplot(aes(x=year, y=purge_freq, fill=region2)) + 
  geom_col( width=0.4) +    # the heights of the bars to represent values in the data
  facet_wrap(~ region2,  nrow = 1)    +
  ylab("Probability of purges per country")+
  xlab("Year") +
  theme(legend.position="none") 


#### Function: Multiplot
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}



### Multiplot (Figure 1 in the paper)

multiplot(p1, p2, cols= 1)  # "trend_rr.pdf"

# high resolution
dev.print(tiff, "C:/Users/Dell User/Dropbox/Military Purge Dataset Project/Paper/trend_rr.tiff", res=600, height=8, width=7, units="in")


# Save directly at 
#"C:/Users/Dell User/Dropbox/Military Purge Dataset Project/Paper/trend_rr.pdf", width=8, height=7)




#####
##### Here I'm Checking the probabilit 1 case in Eastern/Western Europe  

dat %>% filter(gwf_regime!="democracy" & gwf_regime!="foreign/warload")%>%  
  group_by(ccode, year) %>% 
  summarise(
    region2 = first(region2), 
    sum_purge = sum(purge_v1) , 
    leader= first(leader), 
    idacr=first(idacr)
  ) %>% ungroup() %>% 
  mutate(   
    purge = ifelse(sum_purge>0, 1, 0) , 
    region2=ifelse(region2=="Eastern/Western Europe", "Eastern/Western \n Europe", 
                   ifelse(region2=="Latin America", "Latin \n America", 
                          ifelse(region2=="North Africa and the Middle East", "North Africa and \n Middle East", 
                                 ifelse(region2=="Sub-Saharan Africa", "Sub-Saharan \n Africa", region2))))  
    
  ) %>% 
  group_by(region2, year) %>% 
  summarise( 
    a_country = length( unique(ccode)),
    npurge = sum(purge, na.rm=TRUE), 
    purge_freq =npurge/ a_country, 
    fccode=first(idacr), 
    lccode=last(idacr), 
    fleader=first(leader), 
    lleader=last(leader) 
    
  )  %>%    filter(purge_freq==1)  ## checking which year exactly are the 1 probability 
  
#region2                      year a_country npurge purge_freq fccode lccode fleader lleader     
#<chr>                       <dbl>     <int>  <dbl>      <dbl> <chr>  <chr>  <chr>   <chr>       
#  1 "Eastern/Western \n Europe"  2001         2      2          1 RUS    GRG    Putin   Shevardnadze
#2 "Eastern/Western \n Europe"  2004         1      1          1 RUS    RUS    Putin   Putin       
#3 "Eastern/Western \n Europe"  2005         1      1          1 RUS    RUS    Putin   Putin       
 


#### Check GWF 


#GWF book
gwf <- read.dta13(  "C:/Users/pjb13215/Dropbox/DATA/GWF-Autocratic-Regimes-1.2/GWFBookData/GWF.dta") 
gwf <- as_tibble(gwf) 
 
gwf %>% 
    select(cowcode, year, region, gwf_leadername,gwf_country) %>% 
  filter(region=="ceeurope" | region=="weu")  %>% 
  filter(year==2001|year==2004 | year==2005)
 
###  cowcode  year region   gwf_leadername gwf_country
#<int> <int> <chr>    <chr>          <chr>      
# 1     365  2001 ceeurope Putin          Russia     
#2     365  2004 ceeurope Putin          Russia     
#3     365  2005 ceeurope Putin          Russia     
#4     370  2001 ceeurope Lukashenko     Belarus    
#5     370  2004 ceeurope Lukashenko     Belarus    
#6     370  2005 ceeurope Lukashenko     Belarus    
#7     372  2001 ceeurope Shevardnadze   Georgia    

####  2001 and 2004m 2005, Belarus is autocracy, but we did not include in the purge data (it's missing as we couldn't code them)


##### Checking the trends above
## Why the probability and the number of purges look exactly the same for Asia.
dat %>% filter(gwf_regime!="democracy" & gwf_regime!="foreign/warload")%>%  
  group_by(ccode, year) %>% 
  summarise(
    region2 = first(region2), 
    sum_purge = sum(purge_v1) , 
    leader= first(leader)
  ) %>% ungroup() %>% 
  mutate(   
    purge = ifelse(sum_purge>0, 1, 0) , 
    region2=ifelse(region2=="Eastern/Western Europe", "Eastern/Western \n Europe", 
                   ifelse(region2=="Latin America", "Latin \n America", 
                          ifelse(region2=="North Africa and the Middle East", "North Africa and \n Middle East", 
                                 ifelse(region2=="Sub-Saharan Africa", "Sub-Saharan \n Africa", region2))))  
    
  ) %>% 
  group_by(region2, year) %>% 
  summarise( 
    a_country = length( unique(ccode)),
    npurge = sum(purge, na.rm=TRUE), 
    purge_freq =npurge/ a_country, 
    fccode=first(ccode), 
    lccode=last(ccode), 
    fleader=first(leader), 
    lleader=last(leader) 
    
  )  %>%    filter(region2=="Asia")  %>% print(n=Inf) 


###### Calculating the regional average 

dat %>% filter(gwf_regime!="democracy" & gwf_regime!="foreign/warload")%>%  
  group_by(ccode, year) %>% 
  summarise(
    region2 = first(region2), 
    sum_purge = sum(purge_v1) , 
    leader= first(leader)
  ) %>% ungroup() %>% 
  mutate(   
    purge = ifelse(sum_purge>0, 1, 0) , 
    region2=ifelse(region2=="Eastern/Western Europe", "Eastern/Western \n Europe", 
                   ifelse(region2=="Latin America", "Latin \n America", 
                          ifelse(region2=="North Africa and the Middle East", "North Africa and \n Middle East", 
                                 ifelse(region2=="Sub-Saharan Africa", "Sub-Saharan \n Africa", region2))))  
    
  ) %>% 
  group_by(region2, year) %>% 
  summarise( 
    a_country = length( unique(ccode)),
    npurge = sum(purge, na.rm=TRUE), 
    purge_freq =npurge/ a_country, 
    fccode=first(ccode), 
    lccode=last(ccode), 
    fleader=first(leader), 
    lleader=last(leader) 
    
  )  %>% ungroup() %>% 
  group_by(region2)%>% 
  summarise(
    ave_purge = mean(purge_freq)
  )

# region2                           ave_purge
#<chr>                                 <dbl>
#  1 Asia                                  0.251
#2 "Eastern/Western \n Europe"           0.243
#3 "Latin \n America"                    0.188
#4 "North Africa and \n Middle East"     0.250
#5 "Sub-Saharan \n Africa"               0.262






################
################
#### Figure 2 (Main Text)
################
################ 

### leader characteristics and probability of purges 

#### Boxplot (distribution of purge probability)  
### Facet by leaders type (note that i didn't check multiple-leaders year cases yet.)

lead_dat <- dat %>% 
  group_by(obsid) %>%
  summarise(
    count=n(), 
    leaderciv= leaderciv[which(!is.na(leaderciv))[1]],
    leadermil= leadermil[which(!is.na(leadermil))[1]],
    leaderrebel= leaderrebel[which(!is.na(leaderrebel))[1]],  
    all_purge = mean(purge_v1), 
    top_rank= mean(purge_v1_top), 
    violence= mean(violent_purge_v1), 
    topviolence_freq = mean(violentpurge_v1_top),
    gwf_regime= first(gwf_regime),
    leader = first(leader),
    region2=first(region2)
  )  %>% 
  mutate(
    leadertype = ifelse(gwf_regime=="monarchy", "royal", 
                        ifelse(leaderciv==1, "civilian", 
                               ifelse(leadermil==1, "military", 
                                      ifelse(leaderrebel==1, "rebel", NA)))),
    power_length = ifelse(count <15, "Less than 14 yrs", "More than 15 yrs"),
    power_length = factor(power_length, levels = c("Less than 14 yrs", "More than 15 yrs" ))
  )


lead_dat %>% 
  filter(    gwf_regime!="democracy" & gwf_regime!="foreign/warload") %>%  
  count(leadertype)
#civilian     171

#2 military     192

#3 rebel         36

#4 royal         37
# 5 NA             3
## 436 total  
 

#### Reshape the data 

a <-  lead_dat %>% 
  gather(all_purge , top_rank, violence,   key = "purgetype", value="proportion") %>% 
  mutate(purgetype= ifelse(purgetype=="all_purge",  "All purge",
                           ifelse(purgetype=="top_rank", "Top-ranked \npurge", 
                                  ifelse(purgetype=="violence", "Violent \npurge", purgetype)))) %>% 
  filter(   !is.na(leadertype) , count> 1 & gwf_regime!="democracy" & gwf_regime!="foreign/warload") %>%  
  mutate(leadertype=ifelse(leadertype=="royal", "Royal", 
                           ifelse(leadertype=="civilian", "Civilian", 
                                  ifelse(leadertype=="military", "Military",
                                         ifelse(leadertype=="rebel", "Rebel", leadertype))))) %>% 
  ggplot(aes(leadertype, proportion, colour=purgetype)) + 
  geom_boxplot()+ 
  theme(legend.title=element_blank()) +
  theme(legend.position=c(0.8, 0.9))+
  ggtitle("All leaders \n ") +
  theme(plot.title = element_text(size =10)) +
  ylab("Probability of military purge per year") + 
  xlab("Leader type")
 

# Reshape the data 
b <-  lead_dat %>% 
  filter(all_purge !=0) %>%  # only include those leaders who at least have one purge during tenure  
  gather(all_purge, top_rank, violence,     key = "purgetype", value="proportion") %>% 
  mutate(purgetype= ifelse(purgetype=="all_purge",  "All purge",
                           ifelse(purgetype=="top_rank", "Top-ranked \npurge", 
                                  ifelse(purgetype=="violence", "Violent \npurge", purgetype)))) %>% 
  filter(     !is.na(leadertype) , count> 1 & gwf_regime!="democracy" & gwf_regime!="foreign/warload") %>%  
  mutate(leadertype=ifelse(leadertype=="royal", "Royal", 
                           ifelse(leadertype=="civilian", "Civilian", 
                                  ifelse(leadertype=="military", "Military",
                                         ifelse(leadertype=="rebel", "Rebel", leadertype))))) %>% 
  ggplot(aes(leadertype, proportion, colour=purgetype)) + 
  geom_boxplot() +   theme(legend.title=element_blank()) +
  theme(legend.position=c(0.8, 0.9))+
  ggtitle("Leaders with at least one purge \n")+
  theme(plot.title = element_text(size =10)) +
  ylab("Probability of military purge per year") + 
  xlab("Leader type")



c <- lead_dat %>% 
  filter(count >=10) %>%
  gather(all_purge , top_rank, violence,   key = "purgetype", value="proportion") %>% 
  mutate(purgetype= ifelse(purgetype=="all_purge",  "All purge",                                          
                            ifelse(purgetype=="top_rank", "Top-ranked \npurge",                                       
                                    ifelse(purgetype=="violence", "Violent \npurge", purgetype)))) %>% 
  filter(   !is.na(leadertype) , count> 1 & gwf_regime!="democracy" & gwf_regime!="foreign/warload") %>%  
  mutate(leadertype=ifelse(leadertype=="royal", "Royal", 
                           ifelse(leadertype=="civilian", "Civilian", 
                                  ifelse(leadertype=="military", "Military",
                                         ifelse(leadertype=="rebel", "Rebel", leadertype))))) %>% 
  ggplot(aes(leadertype, proportion, colour=purgetype)) + 
  geom_boxplot() +  
  expand_limits(y=c(0,1))+
  scale_y_continuous(breaks=c(0, 0.25, 0.5, 0.75, 1)) +
  theme(legend.title=element_blank()) +
  theme(legend.position=c(0.8, 0.9))+
  ggtitle("Leaders who survived more \nthan 10 years")+
  theme(plot.title = element_text(size =10)) +
  ylab("Probability of military purge per year") + 
  xlab("Leader type")


grid.arrange(a, b, c,  ncol=3)

# save directly "civmil_2rev" (size 7*5,35)


# high resolution
dev.print(tiff, "C:/Users/Dell User/Dropbox/Military Purge Dataset Project/Paper/civmil_2rev.tiff", res=600, height=6, width=8, units="in")








################
############## Figure 3 (in text)
###############


##### The List of Leaders with The Maximum Number of Total Prge 
dat %>%
  group_by(obsid) %>%
  summarise(
    leader =  leader[which(!is.na(leader))[1]],
    idacr = idacr[which(!is.na(idacr))[1]],
    count = n(),
    totalpurge = sum(purge_v1, na.rm=TRUE),
    totalviolence = sum(violent_purge_v1, na.rm=TRUE),
    purge_freq = mean(purge_v1, na.rm=TRUE),
    violence_freq = mean(violent_purge_v1, na.rm=TRUE),
    gwf_regime= first(gwf_regime),
    region2 = first(region2)
  )  %>% 
  mutate(
    leader_ccode = paste(leader, idacr, sep=" ")
    )%>%
  arrange(desc(totalpurge)) %>%
  filter(between(totalpurge, 7, max(totalpurge))) %>%
  filter(gwf_regime!="foreign/warload") %>%
  mutate(
    gwf_regime=ifelse(gwf_regime=="monarchy", "Monarchy", 
                      ifelse(gwf_regime=="party", "Party", 
                             ifelse(gwf_regime=="personal", "Personalist", 
                                    ifelse(gwf_regime=="military", "Military", gwf_regime))) )
  )%>%
  ggplot()+
  geom_bar( stat='identity', aes(x=reorder(leader_ccode, totalpurge),y=totalpurge, fill=gwf_regime ), width=.3) +
  theme(legend.position = "bottom") +  theme(legend.title=element_blank())+
  scale_fill_discrete(breaks=c("Personalist","Military", "Party", "Monarchy" ) ) + 
  ylab("Total number of purges")+
  xlab("") +
  coord_flip()


# high resolution
dev.print(tiff, "C:/Users/Dell User/Dropbox/Military Purge Dataset Project/Paper/maxleader.tiff", res=600, height=7, width=6, units="in")


# Save
## ggsave("C:/Users/Dell User/Dropbox/Military Purge Dataset Project/Paper/maxleader.pdf", width=8, height=5)
# make twice as big as on screen)





 



###############################
###############################
## How to use MPD 
###############################
##############################


library("glmmTMB")
library("bbmle")



## data 
 
dat <- read_rds(   "Purgedata_version2_2019/FINAL/jpr_RR_analysis_v1.rds" ) 

## 3,601 observations 
 

## Exclude "democracy" in GWF (and missing ) - this leads to 3,220 observations 
dat <- dat   %>% 
  filter(gwf_nonautocracy!="democracy") 



## country year

dat %>% filter(!is.na(Government_OneSided_best_fatality_estimate)) %>% 
  count(ccode) %>% dim()   ## 108 

dat %>%  filter(!is.na(Government_OneSided_best_fatality_estimate)) %>% 
  count(year) %>% print(n=Inf)
 


########
######## Model Choide /Model Fit   
########

### Over-dispersion? 
mean(dat$Government_OneSided_best_fatality_estimate)#  171.8314
var(dat$Government_OneSided_best_fatality_estimate) # 77667276
sd(dat$Government_OneSided_best_fatality_estimate) #8812.904
max(dat$Government_OneSided_best_fatality_estimate) # 500000 (5e+05)

# excess zero?
dat %>% count(Government_OneSided_best_fatality_estimate==0) # 25 times more 0 than 1
#FALSE                                               126
# TRUE                                               3094


#### negative binomial 
ng1 <- glmmTMB(Government_OneSided_best_fatality_estimate ~ binary_past5_Lsize_purge_v7 +   
                 supportparty  + log(gwf_leader_duration) + leadermil + leaderciv+  
                 log(lag1_Government_OneSided_best_fatality_estimate +1  )   +   
                 PRIO_civilwar+ 
                 ln_pop_pwt + ln_rgdpech_pwt + 
                 (1|ccode)  ,  data = dat, family=nbinom2
)

summary(ng1 ) 


## Zero-Inflation model with Poisson link 
poisson1 <- update( ng1, . ~ . , 
                    ziformula = ~  lag1_binary_Government_Onesided + binary_past5_Lsize_purge_v7 +  
                      PRIO_civilwar +   ln_pop_pwt +  ln_rgdpech_pwt + (1|ccode) , 
                    family=poisson) 

summary(poisson1) 

## Zero-Inflation model with Negative Binomial Link 

zing1 <- update(poisson1, family=nbinom2)
summary(zing1)


## More control in Zero-Inflation model with Negative Binomial Link   

zing2 <- update(zing1,             
                ziformula= ~ . + supportparty + log(gwf_leader_duration) + leadermil + leaderciv )

summary(zing2)


mstat <- function(model){
  a <- summary(model)
  AIC <- as.numeric(a$AICtab[1])
  BIC <- as.numeric(a$AICtab[2])
  loglike <- as.numeric(a$AICtab[3])
  
  N <- a$nobs 
  return(cbind(N, AIC, BIC, loglike)) 
}

test <- rbind(mstat(ng1), mstat(poisson1),mstat(zing1) ,mstat(zing2) )

test 

AICtab(ng1, poisson1, zing1, zing2)
 
# lower AIC fits better model 
##The cchanges from negative binomial to zero inflation massively improves the model fit. 
## B


 


#####
##### Main Results ( In-Text )
##### 

## binary purge 
m1 <- glmmTMB(Government_OneSided_best_fatality_estimate ~ binary_past5_purge_v7 +   
                supportparty  + log(gwf_leader_duration) + leadermil + leaderciv+  
                log(lag1_Government_OneSided_best_fatality_estimate +1  )   +   
                PRIO_civilwar+ 
                ln_pop_pwt + ln_rgdpech_pwt + 
                (1|ccode) ,  data = dat, family=nbinom2, 
              
              ziformula = ~  lag1_binary_Government_Onesided + binary_past5_purge_v7 +  
                PRIO_civilwar +   ln_pop_pwt +  ln_rgdpech_pwt + (1|ccode) )

summary(m1 ) 

## binary purge 
m2 <- glmmTMB(Government_OneSided_best_fatality_estimate ~ binary_past5_violent_purge_v7 +   
                supportparty  + log(gwf_leader_duration) + leadermil + leaderciv+  
                log(lag1_Government_OneSided_best_fatality_estimate +1  )   +   
                PRIO_civilwar+ 
                ln_pop_pwt + ln_rgdpech_pwt + 
                (1|ccode) ,  data = dat, family=nbinom2, 
              
              ziformula = ~  lag1_binary_Government_Onesided + binary_past5_violent_purge_v7 +  
                PRIO_civilwar +   ln_pop_pwt +  ln_rgdpech_pwt + (1|ccode) )

summary(m2 ) 




##top purge

m3 <- glmmTMB(Government_OneSided_best_fatality_estimate ~ binary_past5_top_purge_v7 +   
                supportparty  + log(gwf_leader_duration) + leadermil + leaderciv+  
                log(lag1_Government_OneSided_best_fatality_estimate +1  )   +   
                PRIO_civilwar+ 
                ln_pop_pwt + ln_rgdpech_pwt + 
                (1|ccode) ,  data = dat, family=nbinom2, 
              
              ziformula = ~  lag1_binary_Government_Onesided + binary_past5_top_purge_v7 +  
                PRIO_civilwar +   ln_pop_pwt +  ln_rgdpech_pwt + (1|ccode) )

summary(m3 ) 


## Large Size (binary)
m4 <- glmmTMB(Government_OneSided_best_fatality_estimate ~ binary_past5_Lsize_purge_v7 +   
                supportparty  + log(gwf_leader_duration) + leadermil + leaderciv+  
                log(lag1_Government_OneSided_best_fatality_estimate +1  )   +   
                PRIO_civilwar+ 
                ln_pop_pwt + ln_rgdpech_pwt + 
                (1|ccode) ,  data = dat, family=nbinom2, 
              
              ziformula = ~  lag1_binary_Government_Onesided + binary_past5_Lsize_purge_v7 +  
                PRIO_civilwar +   ln_pop_pwt +  ln_rgdpech_pwt + (1|ccode) )

summary(m4) 


## Large Size (continuous)
m5 <- glmmTMB(Government_OneSided_best_fatality_estimate ~ past5_Lsize_purge_v7 +   
                supportparty  + log(gwf_leader_duration) + leadermil + leaderciv+  
                log(lag1_Government_OneSided_best_fatality_estimate +1  )   +   
                PRIO_civilwar+ 
                ln_pop_pwt + ln_rgdpech_pwt + 
                (1|ccode) ,  data = dat, family=nbinom2, 
              
              ziformula = ~  lag1_binary_Government_Onesided + binary_past5_Lsize_purge_v7 +  
                PRIO_civilwar +   ln_pop_pwt +  ln_rgdpech_pwt + (1|ccode) )

summary(m5 ) 


AICtab(m1, m2, m3, m4, m5)  


######
######
###### Coefficient Plot 
###### 
######

### FUNCTIONS  
cidat <- function(x, n){  # x sould be the model    
  ci <- confint(x)   # n is the model # 
  
  m <- names(ci[,3 ]) 
  ci <- as_tibble(ci) 
  
  ci$names <- m 
  
  ##ci <- ci %>% filter(str_detect(names, "^cond.")) %>%
  ##  mutate_all(function(x) str_sub(x, 6, str_length(x)))    
  
  ci$model <- factor(n) ## model 1 
  
  ci <- ci %>% 
    rename( lci=    `2.5 %` , 
            hci =   `97.5 %` ) %>% 
    filter(str_detect(names, "^cond.|^zi"))%>% 
    mutate(count = ifelse(str_detect(names, "^cond"), 1, 0), 
           count= factor(count)) %>% 
    mutate(names  = ifelse(count==1, str_sub(names, 6, str_length(names)), 
                           str_sub(names, 4, str_length(names)))
    )
  
  ci$model <- as.character(ci$model) 
  ci$count <- as.character(ci$count)
  
  return(ci) 
}

keycoef <- function(ci){ ## Get the coefficient on the main IV 
  ci <- ci %>% 
    filter(str_detect(names, "past5") &  str_detect(names, "purge")) %>% 
    filter(!str_detect(names, "^zi")) 
  
  return(ci)
} 

#### Multiplot
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}





#######################
### Making Figure (Coefficietn Plots) in Main Text 
#########################

############
### Effects of Differnety Types of Purges (In-Text )
#######
# Apply the above-written functions 
ci1 <- cidat(m1, 1)
ci1<-  keycoef(ci1) 

ci2 <- cidat(m2, 2)
ci2<-  keycoef(ci2) 

ci3 <- cidat(m3, 3)
ci3<-  keycoef(ci3) 

ci4 <- cidat(m4, 4)
ci4<-  keycoef(ci4)   

ci5 <- cidat(m5, 5)
ci5<-  keycoef(ci5)   



# combine all the models 
mdat <- bind_rows(ci1, ci2, ci3, ci4, ci5)   


## Coefficient Plot   
mdat <- 
  mdat %>% mutate(
    names = ifelse(names=="binary_past5_purge_v7", "All purge", 
                   ifelse(names=="binary_past5_violent_purge_v7", "Violent purge",
                          ifelse(names=="binary_past5_top_purge_v7", "Top-ranked purge", 
                                 ifelse(names=="binary_past5_Lsize_purge_v7", "Large-scale purge", 
                                        "Large-scale purge (count)"))))
  ) %>% mutate(
    names = as.factor(names), 
    names = fct_relevel(names, "All purge", "Violent purge", "Top-ranked purge", "Large-scale purge", 
                        "Large-scale purge (count)") , 
    names = fct_rev(names)  ## reverse the order of factor levels 
  )  

levels(mdat$names)  <- gsub(" ", "\n", levels(mdat$names)  )

mdat<-
  mdat %>% mutate(
    model= ifelse(model==1, "Model 1", 
                  ifelse(model==2, "Model 2", 
                         ifelse(model==3, "Model 3", 
                                ifelse(model==4, "Model 4", 
                                       "Model 5")))), 
    model = as.factor(model)
  ) 



## Plot (In-Text Figure)
p1 <- 
  mdat %>% ggplot( aes(names, Estimate, colour=model)) + 
  geom_pointrange(aes(ymin = lci , ymax = hci)) + 
  geom_hline(yintercept = 0, lty=2) +  
  theme_bw()  + 
  coord_flip()   +
  labs(
    title="Effects of military purges",
    x= "Purge variales", 
    y= "", 
    colour= ""
  ) +  
  theme(axis.title.x = element_text(face="bold", colour="#990000", size=20)) +
  scale_colour_discrete(name="Models",
                        breaks=c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5"),
                        labels=c("Model 1 \n ", "Model 2 \n", "Model 3 \n ", "Model 4 \n", "Model 5 \n")) + 
  theme(legend.position="bottom") + 
  theme(legend.title=element_blank()) + 
  guides(col = guide_legend(nrow = 2))




##### 
##### Coefficient Plot for Model 4 (In-Text ) 
##### 
sdat <- cidat(m4, 4)

sdat <- sdat %>%  filter(
  !str_detect(names, "(Intercept)")
)

sdat <- 
  sdat %>% mutate(
    names = ifelse(names=="binary_past5_Lsize_purge_v7", "Large-scale purge", 
                   ifelse(names=="supportparty", "Party", 
                          ifelse(names=="log(gwf_leader_duration)", "Tenure", 
                                 ifelse(names=="leadermil", "Military", 
                                        ifelse(names=="leaderciv", "Civilian", 
                                               ifelse(names=="log(lag1_Government_OneSided_best_fatality_estimate + 1)", "Lag DV", 
                                                      ifelse(names=="PRIO_civilwar", "Civil war", 
                                                             ifelse(names=="ln_pop_pwt", "Population", 
                                                                    ifelse(names=="ln_rgdpech_pwt", "GDP/capita", 
                                                                           ifelse(names=="zi~lag1_binary_Government_Onesided", "Zi:Lag DV", 
                                                                                  ifelse(names=="zi~PRIO_civilwar", "Zi:Civil war", 
                                                                                         ifelse(names=="zi~ln_pop_pwt", "Zi:Population", 
                                                                                                ifelse(names=="zi~binary_past5_Lsize_purge_v7", "Zi:Large scale purge", 
                                                                                                       "Zi:GDP/capita")))))))))))
                   ))
  ) %>% mutate(
    
    names = as.factor(names), 
    names = fct_relevel(names, "Large-scale purge",   "Population", "GDP/capita", "Civil war", 
                        "Lag DV", "Party", "Tenure", "Military", "Civilian" ,  "Zi:Large scale purge", "Zi:Population","Zi:GDP/capita",  "Zi:Civil war",  "Zi:Lag DV"),
    names = fct_rev(names) 
  )


levels(sdat$names)  <- gsub(" ", "\n", levels(sdat$names)  )

sdat <- sdat %>% mutate(
  model = "Model 3", 
  model=as.factor(model)
) %>% mutate(
  count = ifelse(count==1,"Count",  "Zero \n inflation" )
)

## Plot (In-Text Figure)
p2 <- 
  sdat %>% ggplot( aes(names, Estimate, colour=count)) + 
  geom_pointrange(aes(ymin = lci , ymax = hci,   fatten = 1.5)) + 
  geom_hline(yintercept = 0, lty=2) +  
  theme_bw()  + 
  coord_flip()   +
  labs(
    title="Model 4",
    x= "", 
    y= "", 
    colour= ""
  )   + 
  theme(legend.position="bottom")





#### Combine two figure 
multiplot(p1, p2, cols=2)
## save as "coef1.pdf"




# high resolution
dev.print(tiff, "C:/Users/Dell User/Dropbox/Military Purge Dataset Project/Paper/coef1.tiff", res=600, height=7, width=6, units="in")


# Save


 



###########################
###########################
#### Appendix Results : Coefficient Plots 
###########################
###########################


# Apply the above-written functions 
ci1 <- cidat(m1, 1)
ci2 <- cidat(m2, 2)
ci3 <- cidat(m3, 3)
ci4 <- cidat(m4, 4)
ci5 <- cidat(m5, 5)

# combine all the models 
bdat <- bind_rows(ci1, ci2, ci3, ci4, ci5)   


## main models 

bdat <- bdat %>%  filter(
  !str_detect(names, "(Intercept)")
) 

bdat <- bdat %>% mutate(
  names = ifelse(names=="binary_past5_Lsize_purge_v7", "Large-Scale \n Purge", 
                 ifelse(names=="binary_past5_purge_v7" , "All Purge", 
                        ifelse(names=="binary_past5_violent_purge_v7", "Violent Purge", 
                               ifelse(names=="binary_past5_top_purge_v7", "Top-Ranked \n Purge", 
                                      ifelse(names=="past5_Lsize_purge_v7", "Large-Scale \n Purge (Count)", 
                                             ifelse(names=="zi~binary_past5_purge_v7", "Zi:All Purge", 
                                                    ifelse(names=="zi~binary_past5_violent_purge_v7", "Zi:Violent Purge", 
                                                           ifelse(names=="zi~binary_past5_top_purge_v7", "Zi:Top-Ranked \n Purge", 
                                                                  ifelse(names=="supportparty", "Party", 
                                                                         ifelse(names=="log(gwf_leader_duration)", "Tenure", 
                                                                                ifelse(names=="leadermil", "Military", 
                                                                                       ifelse(names=="leaderciv", "Civilian", 
                                                                                              ifelse(names=="log(lag1_Government_OneSided_best_fatality_estimate + 1)", "LagDV", 
                                                                                                     ifelse(names=="PRIO_civilwar", "Civil war", 
                                                                                                            ifelse(names=="ln_pop_pwt", "Population", 
                                                                                                                   ifelse(names=="ln_rgdpech_pwt", "GDP/capita", 
                                                                                                                          ifelse(names=="zi~lag1_binary_Government_Onesided", "Zi:LagDV", 
                                                                                                                                 ifelse(names=="zi~PRIO_civilwar", "Zi:Civil War", 
                                                                                                                                        ifelse(names=="zi~ln_pop_pwt", "Zi:Population", 
                                                                                                                                               ifelse(names=="zi~binary_past5_Lsize_purge_v7", "Zi:Large-Scale \n Purge", 
                                                                                                                                                      "Zi:GDP/capita"))))))))))))
                                                           ))))))))
)

bdat <- bdat %>% mutate(  
  names = as.factor(names), 
  names = fct_relevel(names, "All Purge", "Violent Purge", "Top-Ranked \n Purge",  "Large-Scale \n Purge",  "Large-Scale \n Purge (Count)",   "Population", "GDP/capita", "Civil war", 
                      "LagDV", "Party", "Tenure", "Military", "Civilian" , 
                      "Zi:All Purge", "Zi:Violent Purge", "Zi:Top-Ranked \n Purge", "Zi:Large-Scale \n Purge", 
                      "Zi:Population","Zi:GDP/capita", "Zi:Civil War", "Zi:LagDV"),
  names = fct_rev(names) 
)


#levels(bdat$names)  <- gsub(" ", "\n", levels(bdat$names)  )

bdat <- bdat %>% mutate(
  model = ifelse(model==1, "Model 1", 
                 ifelse(model==2, "Model 2", 
                        ifelse(model==3, "Model 3", 
                               ifelse(model==4, "Model 4", 
                                      "Model 5")))),
  model=as.factor(model), 
  model = fct_relevel(model, "Model 1", "Model 2", "Model 3", "Model 4", "Model 5"), 
  model = fct_rev(model)
) 




####
#### Plot (Appendix Figure 1)
#### 

## Appendix 1 (main result tables )
bdat %>% ggplot( aes(names, Estimate ) ) +  
  geom_pointrange(aes(ymin = lci , ymax = hci,  colour=model, fatten = 1.5), 
                  lwd = 1/2, position = position_dodge(width = 1/2), shape =21 ) + 
  geom_hline(yintercept = 0, lty=2) +  
  theme_bw()  + 
  coord_flip()   + 
  labs(
    title="",
    x= "", 
    y= "", 
    colour= ""
  )   +   
  scale_colour_discrete(name="",
                        breaks=c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5"),
                        labels=c("Model 1 \n (N=2413)", "Model 2  \n (N=2413)", "Model 3 \n (N=2413)", 
                                 "Model 4 \n (N=2413)", "Model 5 \n (N=2413)")
  )+ 
  theme(legend.position="bottom")


## Save as "appcoef1. pdf" 





##############
##############
##### Appendix (Results with more Control variables)
##############
############## 

## Binary Large Size (Model 3 in the main text) 

f1 <- glmmTMB(Government_OneSided_best_fatality_estimate ~ binary_past5_Lsize_purge_v7 +   
                supportparty  + log(gwf_leader_duration) + leadermil + leaderciv+  
                log(lag1_Government_OneSided_best_fatality_estimate +1  )   +   
                PRIO_civilwar+ 
                ln_pop_pwt + ln_rgdpech_pwt + 
                (1|ccode) ,  data = dat, family=nbinom2, 
              
              ziformula = ~  lag1_binary_Government_Onesided + binary_past5_Lsize_purge_v7 +  
                PRIO_civilwar +   ln_pop_pwt +  ln_rgdpech_pwt + (1|ccode) )

summary(f1) 



## With ICC Control: 

f2 <- update(f1,
             . ~ . + ICC_post_ratification )

summary(f2) 

## Aid control -->> No convergence

## Rebel civilian killing (only binary converges): The results same
f3 <- update(f1,
             . ~ . + binary_rebel_onesided) 
summary(f3)

## Ethnic excluded population: the same
f4 <- update(f1,
             . ~ . + EPR_exclpop) 
summary(f4)


## Protest (all good, but need to separate violent and nonviolent to be converged)
f5 <- update(f1, 
             . ~ . + navco_violent + navco_nonviolent)

summary(f5)


## personalism
f6 <- update(f1,
             . ~ . + latent_personalism)

summary(f6)


## International conflict 
f9 <- update(f1,
             . ~ . + PRIO_interstate_conflict)

summary(f9)



## Year-random effects in addtition to Country-random effects (results same)
f7 <-   glmmTMB(Government_OneSided_best_fatality_estimate ~ binary_past5_Lsize_purge_v7 +   
                  supportparty  + log(gwf_leader_duration) + leadermil + leaderciv+  
                  log(lag1_Government_OneSided_best_fatality_estimate +1  )   +   
                  PRIO_civilwar+ 
                  ln_pop_pwt + ln_rgdpech_pwt + 
                  (1|year) + (1|ccode) ,  data = dat, family=nbinom2, 
                
                ziformula = ~  lag1_binary_Government_Onesided + binary_past5_Lsize_purge_v7 +  
                  PRIO_civilwar +   ln_pop_pwt +  ln_rgdpech_pwt + (1|year) + (1|ccode) )

summary(f7) 





#######################
## FUNCTIONS :::::
#######################

redat <- function(dat){
  
  dat <- dat %>%  filter(
    !str_detect(names, "(Intercept)")
  ) 
  
  dat <- dat %>% mutate(
    names = ifelse(names=="binary_past5_Lsize_purge_v7", "Large-Scale \n Purge", 
                   ifelse(names=="binary_past5_purge_v7" , "All Purge", 
                          ifelse(names=="binary_past5_violent_purge_v7", "Violent Purge", 
                                 ifelse(names=="past5_Lsize_purge_v7", "Large-Scale \n Purge (Count)", 
                                        ifelse(names=="zi~binary_past5_purge_v7", "Zi:All Purge", 
                                               ifelse(names=="zi~binary_past5_violent_purge_v7", "Zi:Violent Purge", 
                                                      ifelse(names=="supportparty", "Party", 
                                                             ifelse(names=="log(gwf_leader_duration)", "Tenure", 
                                                                    ifelse(names=="leadermil", "Military", 
                                                                           ifelse(names=="leaderciv", "Civilian", 
                                                                                  ifelse(names=="log(lag1_Government_OneSided_best_fatality_estimate + 1)", "LagDV", 
                                                                                         ifelse(names=="PRIO_civilwar", "Civil war", 
                                                                                                ifelse(names=="ln_pop_pwt", "Population", 
                                                                                                       ifelse(names=="ln_rgdpech_pwt", "GDP/capita", 
                                                                                                              ifelse(names=="zi~lag1_binary_Government_Onesided", "Zi:LagDV", 
                                                                                                                     ifelse(names=="zi~PRIO_civilwar", "Zi:Civil War", 
                                                                                                                            ifelse(names=="zi~ln_pop_pwt", "Zi:Population", 
                                                                                                                                   ifelse(names=="zi~binary_past5_Lsize_purge_v7", "Zi:Large-Scale \n Purge",
                                                                                                                                          ifelse(names=="zi~ln_rgdpech_pwt", "Zi:GDP/capita", 
                                                                                                                                                 names
                                                                                                                                          )))))))))))
                                                             ))))))))
  ) 
  
  return(dat)
} 



# Apply the above-written functions 
ci1 <- cidat(f1, 1)
ci2 <- cidat(f2, 2)
ci3 <- cidat(f3, 3)
ci4 <- cidat(f4, 4)
ci5 <- cidat(f5, 5) 
ci6 <- cidat(f6, 6) 

# combine all the models 
bdat <- bind_rows(ci1, ci2, ci3, ci4, ci5, ci6)   

## apply function 
bdat <- redat(bdat)

### adjust new control variables 
advar <- function(bdat){
  bdat <- bdat %>% mutate(
    names=ifelse(names=="binary_rebel_onesided", "Rebel civilian \n killing", 
                 ifelse(names=="ICC_post_ratification", "ICC \n Ratification", 
                        ifelse(names=="EPR_exclpop", "Ethnic \n exclusion", 
                               ifelse(names=="navco_violent", "Violent \n protest", 
                                      ifelse(names=="navco_nonviolent", "Nonviolent \n protest", 
                                             ifelse(names=="latent_personalism", "Personalist", 
                                                    names))))))
  )
  return(bdat)
}

bdat <- advar(bdat)

bdat <- bdat %>% mutate(  
  names = as.factor(names), 
  names = fct_relevel(names,   "Large-Scale \n Purge",   "Population", "GDP/capita", "Civil war", 
                      "LagDV", "Party", "Tenure", "Military", "Civilian" , 
                      "ICC \n Ratification",  "Rebel civilian \n killing", "Ethnic \n exclusion", 
                      "Violent \n protest",  "Nonviolent \n protest",  "Personalist", 
                      "Zi:Large-Scale \n Purge", 
                      "Zi:Population","Zi:GDP/capita", "Zi:Civil War", "Zi:LagDV"),
  names = fct_rev(names) 
)


#levels(bdat$names)  <- gsub(" ", "\n", levels(bdat$names)  )

bdat <- bdat %>% mutate(
  model = ifelse(model==1, "Model 1", 
                 ifelse(model==2, "Model 2", 
                        ifelse(model==3, "Model 3", 
                               ifelse(model==4, "Model 4", 
                                      ifelse(model==5, "Model 5", 
                                             "Model 6" ))))), 
  model=as.factor(model), 
  model = fct_relevel(model, "Model 1", "Model 2", "Model 3", "Model 4", "Model 5", "Model 6"), 
  model = fct_rev(model)
) 




####
#### Plot (Appendix Figure 1)
#### 

## Appendix 1 (main result tables )
bdat %>% ggplot( aes(names, Estimate ) ) +  
  geom_pointrange(aes(ymin = lci , ymax = hci,  colour=model, fatten = 1.5), 
                  lwd = 1/2, position = position_dodge(width = 1/2), shape =21 ) + 
  geom_hline(yintercept = 0, lty=2) +  
  theme_bw()  + 
  coord_flip()   + 
  labs(
    title="",
    x= "", 
    y= "", 
    colour= ""
  )   +   
  scale_colour_discrete(name="",
                        breaks=c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5", "Model 6"),
                        labels=c("Model 1 \n (N=2413)", "Model 2  \n (N=2413)", "Model 3 \n (N=2413)", "Model 4 \n (N=2368)",
                                 "Model 5 \n (N=2413)", "Model 6  \n (N=2413)"  )
  )+ 
  theme(legend.position="bottom")




## Save (appcoef2.pdf)




##############
##############
##### Appendix 3 (Count of Purge IV and Control variables)
##############
############## 

## Count Large Size (Model 4 in the main text) 

f1 <- glmmTMB(Government_OneSided_best_fatality_estimate ~  past5_Lsize_purge_v7 +   
                supportparty  + log(gwf_leader_duration) + leadermil + leaderciv+  
                log(lag1_Government_OneSided_best_fatality_estimate +1  )   +   
                PRIO_civilwar+ 
                ln_pop_pwt + ln_rgdpech_pwt + 
                (1|ccode) ,  data = dat, family=nbinom2, 
              
              ziformula = ~  lag1_binary_Government_Onesided + binary_past5_Lsize_purge_v7 +  
                PRIO_civilwar +   ln_pop_pwt +  ln_rgdpech_pwt + (1|ccode) )

summary(f1) 



## With ICC Control: 

f2 <- update(f1,
             . ~ . + ICC_post_ratification )

summary(f2) 

## Aid control -->> No convergence

## Rebel civilian killing (only binary converges): The results same
f3 <- update(f1,
             . ~ . + binary_rebel_onesided) 
summary(f3)

## Ethnic excluded population: the same
f4 <- update(f1,
             . ~ . + EPR_exclpop) 
summary(f4)


## Protest (all good, but need to separate violent and nonviolent to be converged)
f5 <- update(f1, 
             . ~ . + navco_violent + navco_nonviolent)

summary(f5)


## personalism
f6 <- update(f1,
             . ~ . + latent_personalism)

summary(f6)



# Apply the above-written functions 
ci1 <- cidat(f1, 1)
ci2 <- cidat(f2, 2)
ci3 <- cidat(f3, 3)
ci4 <- cidat(f4, 4)
ci5 <- cidat(f5, 5) 
ci6 <- cidat(f6, 6) 

# combine all the models 
bdat <- bind_rows(ci1, ci2, ci3, ci4, ci5, ci6)   

## apply function 
bdat <- redat(bdat)

bdat <- advar(bdat) ## using the above function 


bdat <- bdat %>% mutate(  
  names = as.factor(names), 
  names = fct_relevel(names,   "Large-Scale \n Purge (Count)",   "Population", "GDP/capita", "Civil war", 
                      "LagDV", "Party", "Tenure", "Military", "Civilian" , 
                      "ICC \n Ratification",  "Rebel civilian \n killing", "Ethnic \n exclusion", 
                      "Violent \n protest",  "Nonviolent \n protest",  "Personalist", 
                      "Zi:Large-Scale \n Purge", 
                      "Zi:Population","Zi:GDP/capita", "Zi:Civil War", "Zi:LagDV"),
  names = fct_rev(names) 
)


#levels(bdat$names)  <- gsub(" ", "\n", levels(bdat$names)  )

bdat <- bdat %>% mutate(
  model = ifelse(model==1, "Model 1", 
                 ifelse(model==2, "Model 2", 
                        ifelse(model==3, "Model 3", 
                               ifelse(model==4, "Model 4", 
                                      ifelse(model==5, "Model 5", 
                                             "Model 6" ))))), 
  model=as.factor(model), 
  model = fct_relevel(model, "Model 1", "Model 2", "Model 3", "Model 4", "Model 5", "Model 6"), 
  model = fct_rev(model)
) 




####
#### Plot (Appendix Figure 1)
#### 

## Appendix 1 (main result tables )
bdat %>% ggplot( aes(names, Estimate ) ) +  
  geom_pointrange(aes(ymin = lci , ymax = hci,  colour=model, fatten = 1.5), 
                  lwd = 1/2, position = position_dodge(width = 1/2), shape =21 ) + 
  geom_hline(yintercept = 0, lty=2) +  
  theme_bw()  + 
  coord_flip()   + 
  labs(
    title="",
    x= "", 
    y= "", 
    colour= ""
  )   +   
  scale_colour_discrete(name="",
                        breaks=c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5", "Model 6"),
                        labels=c("Model 1 \n (N=2413)", "Model 2  \n (N=2413)", "Model 3 \n (N=2413)", "Model 4 \n (N=2368)",
                                 "Model 5 \n (N=2413)", "Model 6  \n (N=2413)"  )
  )+ 
  theme(legend.position="bottom")





## Save as "appcoef3.pdf 













###
### Predicted Value 
###

## My conclusion at this stage (February 2020) is that there is no established way to computer the SE or Confidence intervals for the 
## predicted value. I think something odd translating the zero-inflation compoenetn into prediction. Therefore I use the fixed effects of the 
##conditional model onlyl below.


## https://strengejacke.github.io/ggeffects/articles/randomeffects.html

library("ggeffects") 

#the fixed effects of the conditional model only (type = "fe")
#the fixed effects and zero-inflation component (type = "fe.zi")
#the fixed effects of the conditional model only (population-level), taking the random-effect variances into account (type = "re")
#the fixed effects and zero-inflation component (population-level), taking the random-effect variances into account (type = "re.zi")


ggpredict(m4, terms="binary_past5_Lsize_purge_v7", type="re", 
          condition = c(lag1_binary_Government_Onesided=0, supportparty  = median(dat$supportparty, na.rm=TRUE) , 
                        gwf_leader_duration= mean(dat$gwf_leader_duration, na.rm=TRUE) , 
                        leadermil=median(dat$leadermil, na.rm=TRUE), 
                        leaderciv= median(dat$leaderciv, na.rm=TRUE), 
                        lag1_Government_OneSided_best_fatality_estimate = mean(dat$lag1_Government_OneSided_best_fatality_estimate, na.rm=TRUE), 
                        PRIO_civilwar = median(dat$PRIO_civilwar, na.rm=TRUE)  , 
                        ln_pop_pwt = mean(dat$ln_pop_pwt, na.rm=TRUE) , 
                        ln_rgdpech_pwt  = mean(dat$ln_rgdpech_pwt, na.rm=TRUE), 
                        lag1_binary_Government_Onesided = median( dat$lag1_binary_Government_Onesided, na.rm=TRUE) 
          ) 
)












#####################
##### End      ######
#####################
#################


